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1.  Introduction 

The  existing  lithium  ion  battery  model  in  COSMOL  Inc.  Multi¬ 
physics  3.5a  is  extended  here  by  adding  an  energy  balance  and  the 
temperature  dependence  of  properties  of  the  battery. 

This  thermal  model  is  developed  based  on  the  pseudo  two- 
dimensional  (P2D)  model  which  was  described  in  [1,2]  and  a 
thermal,  electrochemistry  coupled  model.  The  diffusion  coefficient 
of  Li  ions  in  the  solid  phase  and  electrolyte,  the  reaction  rate  con¬ 
stants  of  the  electrochemical  reactions,  the  open  circuit  potentials 
and  the  thermal  conductivity  of  the  binary  electrolyte  depend  on 
the  temperature  in  the  model  presented  here. 

2.  Mathematical  model 


where  i=p  for  the  positive  electrode  and  i  =  n  for  the  negative  elec¬ 
trode.  At  the  center  of  the  particle,  there  is  no  flux,  that  is: 

dcs  i 

at  r  =  0  :  -Dsi  =  0 
5,1  dr 

On  the  surface  of  the  particle,  the  flux  is  equal  to  the  consum¬ 
ing/producing  rate  of  Li  ions  due  to  the  electrochemical  reaction 
occurring  at  the  solid/liquid  interface: 

dcs  i 

at  r  =  Rs j  :  —  Ds  j  —  =  Jj 

where  J;  is  the  flux  of  lithium  ions  away  from  the  surface  of  the 
spherical  particles. 

The  material  balance  for  the  binary  electrolyte  in  the  liquid 
phase  is  given  by: 


A  schematic  of  a  lithium  ion  battery  is  shown  in  Fig.  1.  Gen¬ 
erally,  a  lithium  ion  battery  consists  of  the  current  collector,  the 
positive  electrode,  the  separator  and  the  negative  electrode.  A  lithi- 
ated  organic  solution  fills  the  porous  components  and  serves  as  the 
electrolyte. 

The  material  balance  for  the  Li  ions  in  an  active  solid  material 
particle  is  governed  by  Fields  second  law  in  spherical  coordinates: 


dcs,i  _  n  2  d 

~dF  ~  dr 


(i) 
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d_ 

dx 


DefU  dx 


+  (1  -t°+)aJi 


(2) 


where  i=p ,  s,  and  n  and  a2  is  the  electrode  surface  area  per  unit 
volume  of  the  electrode.  In  the  separator  the  pore  wall  flux  Js  is 
equal  to  zero.  At  the  two  ends  of  the  cell  in  the  x-direction,  there  is 
no  mass  flux,  that  is: 


-D 

-D 


dcp 
eff,p~9x 
dcn 
eff’n  dx 


=  0 

x=0 


x=Lp+Ls+Ln 


=  0 


At  the  interfaces  between  the  positive  electrode/separator 
and  separator/negative  electrode,  the  concentration  of  the  binary 
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Nomenclature 

as 

specific  interfacial  area,  m2  m-3 

c 

concentration  of  the  binary  electrolyte,  mol  m-3 

Cp 

specific  heat,  J  kg-1  K-1 

Cs 

concentration  of  lithium  ion  in  solid,  molm-3 

De 

salt  diffusion  coefficient,  m2  s_1 

Ds 

diffusion  coefficient  of  lithium  ion  in  solid  electrode 
particles,  m2  s-1 

FceW 

cell  voltage,  V 

F 

Faraday’s  constant,  96,487  Ceq-1 

h 

effective  heat  transfer  coefficient,  W  m-2  K-1 

Japp 

applied  current  density,  Am-2 

J 

pore  wall  flux  of  lithium  ions,  molm-2  s_1 

k 

electrochemical  reaction  rate  constant, 

m2-5  mol-0-5  s_1 

L 

thickness  of  battery  component,  m 

R 

gas  constant,  8.3145  J  mol-1 1<-1 

Rs 

radius  of  electrode  particle,  m 

t 

time,  s 

T 

temperature,  K 

t+ 

transference  number  of  lithium  ion 

U 

open-circuit  potential,  V 

Ui,  ref 

open  circuit  potential  of  electrode  i  under  the  refer¬ 
ence  temperature,  V 

X 

spatial  coordinate 

£ 

porosity 

£f 

volume  fraction  of  fillers 

P 

density,  kg  m-3 

electronic  conductivity  of  solid  matrix,  S  m-1 

A 

thermal  conductivity,  Wm-1  K-1 

01 

potential  in  the  solid  phase,  V 

02 

potential  in  the  electrolyte,  V 

Subscripts 

eff 

effective 

max 

maximum 

n 

negative  electrode 

P 

positive  electrode 

s 

separator 

electrolyte  and  its  flux  are  continuous 


-D 

eff-p  3x 

-D  — 
Deff  s  dx 


x=L~ 


~  °effs3x 


=  -D, 


eff,n 


9cn 


x=(Lp+Ls ) 


9x 


x=(Lp+Ls)+ 


The  effective  diffusion  coefficient  of  the  binary  electrolyte  in  the 
liquid  phase  is  corrected  by  porosity: 

n  „  .  _  n  .c  bruggj 


where  the  diffusion  coefficient  of  the  binary  electrolyte  (De  l)  is 
given  in  Ref.  [3]  as  follows: 

De  i  =  1(T4  x  1  o — 4*43 — (54/C7^ — 229 — 5.°x  10_ 3c:x)) — o.22x l0_ 3cz-  ^ 


The  charge  balance  in  the  solid  phase  is  governed  by  Ohm’s  law: 


^eff ,  i 


&4>u 

dx2 


=  aM 


(4) 


where  i=p  and  n.  The  effective  conductivity  is  determined 
by  o eff  j  =  cTj(  1  —  £j  —  E/i ).  The  specific  area  can  be  calculated  by 
a2  =  (3/RS  j)(l  -  £[  -  Sff).  At  the  interface  of  the  current  collector  and 
the  positive  electrode,  the  charge  flux  is  equal  to  the  current  density 
applied  to  the  cell: 


&eff,p 


901, p 
dx 


~  W 

x=0 


where  /app  is  the  applied  current  density. 

At  the  interfaces  of  the  positive  electrode/separator  and  the  sep¬ 
arator/negative  electrode,  there  is  no  charge  flux: 


90i,  p 
CTeff’p  3x 

9<Ai,n 

CTeff’n  dx 


=  0 


x=Lp 


=  0 


x=Lp  +LS 


The  potential  of  the  solid  phase  at  the  right  end  of  the  cell  is  set 
to  zero  <p\,n\x=Lp+Ls+Ln  =  an(^  t^e  potential  of  the  solid  phase  at 
x  =  0,  0i,p|x=o  is  equal  to  £cen  (the  cell  voltage). 

The  charge  balance  in  the  liquid  phase  is  based  on  Ohm’s  law 
and  is  given  by: 


d(  302, ,-V  2RT(l-t°)  3 

dx  yeff9  dx  )+  F  dx 


9(lnq)\ 
3x  J 


=  aith 


(5) 


Cp\x=L~  ~  Cs\x=L+ 
Cs\x=(Lp+Lsr  =C"I  x=(Lp+Ls)+ 


X 

— ► 


Positive  Separator  Negative 
Electrode  Electrode 

Fig.  1.  Schematic  of  a  lithium  ion  battery. 


where  i=p,s  and  n  and  the  specific  conductivity  of  the  binary  elec¬ 
trolyte  is  a  function  of  temperature  and  the  concentration  of  the 
binary  electrolyte  in  the  liquid  phase: 

Keff,i=Ki£ibmSSi 

where  the  expression  of  the  ionic  conductivity  for  the  binary  elec¬ 
trolyte  is  given  in  Ref.  [3]  as  follows: 

x,  =  10“4  x  c,(-10.5  +  0.668  x  10“3c,  +  0.494  x  W~6cf 

+  0.074T— 1.78  x  10_5c,T  —  8.86  x  10“10C?T  -  6.96  x  10“5T2 

+2.80  x  10“8c,T2)2  (6) 


At  the  two  ends  of  the  cell,  there  is  no  charge  flux  in  the  liquid 
phase: 


_K 

eff’P  dx 


=  0  and  -/ceff  n 


902,  n 


x=0 


dx 


=  0 


X=Lp+Ls+Ln 


The  potential  in  the  solution  phase  and  its  flux  are  continuous 
at  the  interfaces  of  the  electrodes  and  the  separator. 
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In  the  above  equations  the  pore  wall  flux,  Jit  is  determined  by 
the  Bulter-Volmer  equation: 

j  /,  ( r  _  r  \0.5  0.5  r0.5 

Ji  —  K-i\Ls,i,  max  Ls,i,  surf;  Ls,i, surrT 


/0.5  F  \ 

(  °-5F  M 

X 

[exp  VkT^) 

and  the  overpotential  of  the  intercalation  reaction  is  given  by: 

*7 i  =  01, i  ~  02, i  ~  Ui 

The  energy  balance  is  given  by  [4,5]: 

dr  02  t 

pCp  ^  0^2"  Qrxn  +  Qrev  +  Qohm  (8) 

with  the  boundary  conditions  determined  by  Newton’s  cooling  law: 


=  h(TOQ-T ) 

x=0 

=  h(T-TOQ) 

x=Lp+Ls+Ln 

where  h  is  the  heat  transfer  coefficient,  is  the  environmental 
temperature,  Qrxn  is  the  total  reaction  heat  generation  rate,  Qrev  is 
the  total  reversible  heat  generation  rate,  Qohm  is  the  total  ohmic 
heat  generation  rate.  The  heat  fluxes  are  defined  by: 

Qrxn  =  Fa/(01  —  02  —  IT) 

M  c  ^dU 
Qrev  —  FaJT  — — 


,  dT 
Xdx 
,  dr 
Xdx 


Qohm  —  ^eff 


(  d<p2\  2/ceffRT  0  d(ln  c)  d(p2 

+Keff(^rJ  +~ — ar 


The  temperature  dependent  open  circuit  potential  of  electrode 
i  is  approximated  by  Taylor’s  first  order  expansion  around  a  refer¬ 
ence  temperature: 

Ui  =  U1.ref  +  (T-Tref)[^]j 

where  U{  ref  is  the  open  circuit  potential  under  the  reference  tem¬ 
perature  Tref. 


3.  The  COMSOL  Multiphysics  (MP)  3.5a  model 

The  mathematical  model  described  in  Section  2  is  a  multi-scale 
model.  We  developed  several  geometries  using  this  software:  a 
ID  geometry  which  consists  of  three  sequentially  connected  lines 
to  represent  the  positive  electrode,  the  separator  and  the  neg¬ 
ative  electrode,  respectively,  a  2D  geometry  which  consists  of 
two  rectangles  to  denote  the  solid  phase  in  the  electrodes.  The 
geometries  are  shown  in  Fig.  2.  The  vertical  coordinate  in  the  2D 
geometry  indicates  the  radial  direction  of  the  solid  particles.  Since 
we  ignore  the  diffusion  of  Li  ions  in  the  x-direction  in  the  parti¬ 
cle,  the  corresponding  diffusion  coefficient  is  set  to  zero  in  this 
direction.  The  concentration  of  the  binary  electrolyte,  the  poten¬ 
tial  in  the  electrolyte,  the  potential  in  the  solid  phase  and  the 
pore  wall  flux  are  solved  in  the  ID  geometry.  The  concentration 
of  Li  ions  in  the  solid  phase  is  solved  in  the  2D  geometry.  The 
pore  wall  flux  is  extruded  from  the  ID  domain  and  projected 
to  the  top  boundary  of  the  2D  geometry  by  using  “subdomain 
extrusion  coupling  variables”  in  COMSOL  Multiphysics.  The  con¬ 
centration  of  Li  ions  on  the  top  boundary  in  the  2D  geometry  is 
projected  to  the  ID  domain  by  using  “boundary  extrusion  coupling 
variables”. 


Fig.  2.  Geometries  and  variables  coupling  between  the  geometries  using  COSMOL 
Multiphysics  3.5a.  (The  top  ID  geometry  consists  of  three  segments  which  denote 
the  positive  electrode,  the  separator  and  the  negative  electrode.  The  bottom  two 
rectangles  represent  the  solid  phases  in  the  positive  electrode  and  the  negative 
electrode,  respectively.  The  pore  wall  flux  is  extracted  from  the  1 D  geometry  and 
projected  to  the  top  boundary  of  the  2D  geometry.  The  concentration  of  Li  ions  on 
the  top  boundary  of  the  2D  geometry  is  projected  to  the  ID  domain  as  the  surface 
concentration  of  the  solid  particles.) 


The  thermal  behavior  of  the  Li  ion  battery  during  pulse  discharge 
is  also  simulated  in  COSMOL  Multiphysics.  The  battery  is  discharged 
for  3000  s  at  C/2  rate  first  and  then  discharged  at  3  C  rate  until  the 
cell  voltage  drops  to  2.5  V.  The  change  of  the  applied  current  density 
is  implemented  by  using  the  smoothed  Heaviside  function  “flsmhs” 
and  is  shown  in  Fig.  3.  The  model  parameters  can  be  found  in  Ref. 
[5].Simulation  results 

Fig.  4  shows  the  temperature  on  the  cell  surface  at  1  C  discharge 
process  under  three  different  cooling  conditions  where  the  heat 
transfer  coefficient  is  1 0.0, 1 .0  and  0.1  W  m-2 1<-1 ,  respectively,  and 
two  limiting  conditions:  the  isothermal  condition  and  the  adiabatic 
condition. 

The  thermal  effect  on  the  cell  voltage  shown  in  Fig.  5.  The  cell 
provides  more  discharge  capacity  when  it  is  placed  in  a  better  heat 
isolation  environment  (i.e.  adiabatic  condition).  In  a  better  isolated 
environment,  the  cell  temperautre  increases  faster  during  the  1  C 
discharge  process  which  results  in  the  higher  diffusion  coefficient 
for  the  binary  electrolyte  and  reduces  the  diffusion  limitations. 

The  reduction  of  the  diffusion  limit  in  the  binary  electrolyte  can 
be  verified  by  comparing  the  concentration  profile  of  the  electrolyte 
under  different  cooling  conditions.  Fig.  6  shows  the  concentration 
profiles  of  the  electrolyte  at  the  end  of  1  C  discharge  process  under 
the  two  limiting  conditions.  The  concentration  profile  under  the 
adiabatic  condition  is  flatter  than  that  in  the  isothermal  case,  which 
indicates  a  better  diffusion  property  in  the  electrolyte  under  the 
adiabatic  condition  than  under  the  isothermal  condition. 
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Fig.  4.  Temperature  on  the  cell  surface  during  1  C  discharge  process  under  different  cooling  conditions. 


Fig.  5.  Cell  voltage  for  1  C  discharge  process  under  different  cooling  conditions. 


Fig.  7  shows  the  cell  temperature  during  the  1  C  discharge  pro¬ 
cess  at  different  current  rates  as  the  heat  transfer  coefficient  is 
1.0Wm_2K_1.  As  expected,  the  cell  gets  hotter  as  the  dishcarge 
current  rate  increases.  It  is  also  noticed  that  the  wave  part  which 
appears  in  beginning  of  the  temperature  curve  at  low  current 
rate  (less  than  2C)  does  not  exist  in  the  high  current  rate  cases. 
The  wave  part  on  the  temperature  curve  is  characterised  by  the 
reversible  heat  generation  during  discharging.  Under  low  current 


Fig.  6.  Comparison  of  the  concentration  profiles  of  the  binary  electrolyte  at  the 
end  of  the  1  C  discharge  process  under  the  isothermal  condition  and  the  adiabatic 
condition. 


rate  discharging,  the  reversible  heat  is  roughly  equivelant  to  the 
ohmic  heat,  but  becomes  unimportant  as  the  discharge  current  rate 
increases. 

The  P2D  model  mentioned  in  Section  2  is  also  useful  for  simulat¬ 
ing  the  discharge  process  with  pulse.  Fig.  8  shows  the  cell  voltage 
during  the  C/2  discharge  for  3000  s  followed  by  a  3  C  pulse  discharge 
until  the  cell  voltage  drops  to  2.5  V.  The  corresponding  temperature 


Fig.  7.  Temperature  on  the  cell  surface  during  discharge  process  under  different 
current  rates  while  the  heat  transfer  coefficient  h  =  1 .0  W  m-2  K-1  where  the  DOD  is 
defined  as:  DOD  =  time  *  C  rate/3600. 
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on  the  surface  of  the  cell  is  also  plotted  in  Fig.  9.  The  surface  tem¬ 
perature  at  the  end  of  the  3  C  pulse  is  slightly  less  than  that  in  the 
pure  3  C  discharge  process. 

Fig.  10  shows  the  concentration  of  the  binary  electrolyte  at 
the  two  ends  of  the  cell  during  the  pulse  discharge  process.  At 


Fig.  10.  Concentration  of  the  binary  electrolyte  at  the  two  ends  of  the  cell  (the 
bottom  line  is  at  x  =  0,  the  top  line  is  at  x  =  1 .65  x  1 0-4  m)  in  the  C/2  discharge  process 
with  a  3  C  pulse.  Mathematical  modeling  of  a  lithium  ion  battery  with  thermal  effects 
in  COMSOL. 


the  beginning  of  the  pulse,  the  concentration  of  the  electrolyte 
changes  extremely,  after  that  it  relaxes  and  tend  to  a  stable 
value. 

5.  Conclusions 

The  thermal  behavior  of  a  lithium  ion  battery  during  galvano- 
static  discharge  w/o  pulse  can  be  predicted  by  using  COMSOL 
Multiphysics  (MP)  3.5a.  Though  better  thermal  isolation  environ¬ 
ment  improves  the  discharge  capacity,  the  higher  cell  temperature 
raises  the  risk  of  thermal  runaway  and  more  rapid  degradation  of 
the  cell. 
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